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In order to perform accurate and stable long-time numerical integration of the Einstein equation, 
several hyperbolic systems have been proposed. We here present numerical comparisons between 
weakly hyperbolic, strongly hyperbolic, and symmetric hyperbolic systems based on Ashtekar's con- 
nection variables. The primary advantage for using this connection formulation in this experiment 
, is that we can keep using the same dynamical variables for all levels of hyperbolicity. Our numerical 

■ code demonstrates gravitational wave propagation in plane symmetric spacetimes, and we compare 

the accuracy of the simulation by monitoring the violation of the constraints. By comparing with 
results obtained from the weakly hyperbolic system, we observe the strongly and symmetric hyper- 
CIh' bolic system show better numerical performance (yield less constraint violation), but not so much 

U ' difference between the latter two. Rather, we find that the symmetric hyperbolic system is not 

always the best in numerical performances. 
\^ ' This study is the premier to present full numerical simulations using Ashtekar's variables. We 

also describe our procedures in detail. 

(N 

^ ' gr-qc/0005003, revised version 

cn ■ 

o : 
o ■ 

^ ■ I. INTRODUCTION 

o : 
o . 

. Numerical relativity - solving the Einstein equation numerically - is now an essential field in gravity research. 

' As is well known, critical collapse in gravity systems was first discovered by numerical simulation Q. The current 
O^" mainstream of numerical relativity is to demonstrate the final phase of compact binary objects related to gravitational 

. wave observations [|| , and these efi^orts are now again shedding light on the mathematical structure of the Einstein 
^P' equations. 

J> . Up to a couple of years ago, the standard Arnowitt-Deser-Misner (ADM) decomposition of the Einstein equation 
, was taken as the standard formulation for numerical relativists. Difficulties in accurate/stable long-term evolutions 
^ ■ were supposed to be overcome by choosing proper gauge conditions and boundary conditions. Recently, however, 
. 5^ 1 several numerical experiments show that the standard ADM is not the best formulation for numerics, and finding a 
better formulation has become one of the main research topics. Q 

One direction in the community is to apply conformally decoupled and tracefree re-formulation of ADM system 
which were first used by Nakamura et al. js) . The usefulness of this re- formulation were confirmed by another groups 
to show a long-term stable numerical evolutions. |^,^. Although there is an effort to show why this re-formulation is 
better than ADM Q, we do not yet know this method is robust for all situations. 

Another alternative approach to ADM is to formulate the Einstein equations to reveal hyperbolicity . A certain 
kind of hyperbolicity of the dynamical equations is essential to analyze their propagation features mathematically, 
and are known to be useful in numerical approximations (we explain these points in §2). The propagation of the 
original ADM constraint equations obeys well-posed behavior but the dynamical equations of the ADM system 
are not a hyperbolic system at all (and these facts can be also applied to the conformally decoupled version) . Several 



^ Note that we are only concerned with the free evolution system of the initial data, that is, we only solve the constraint 
equations on the initial hypersurface. The accuracy and/or stability of the system is normally observed by monitoring the 
violation of constraints during the free evolution. 
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hyperbolic formulations have been proposed to re-express the Einstein equation, with different levels: weakly, strongly 
and symmetric hyperbolic systems (we will discuss in detail in §2). Several numerical tests were also performed in 
this direction, and we can see advantages in numerical stability over the original ADM system (e.g. tests Q of 
Bona-Masso's symmetrizable form tests of Choquet-Bruhat and York (95) 's symmetrizable form but 
the appearance of coordinate shocks is also reported ( in the system of Q ) . A symmetric hyperbolic system of 
, on the other hand, is numerically studied in the context of "conformal Einstein" approach . 

The following questions, therefore, naturally present themselves (cf. Q): (1) Does hyperbolicity actually contribute 
to the numerical accuracy/stability? (2) If so, which level of hyperbolic formulation is practically useful for numerical 
applications? (or does the symmetric hyperbolicity solve all the difficulties?) (3) Are there any other approaches to 
improve the accuracy/stability of the system? 

In this paper, we try to answer these questions with our simple numerical experiments. Such comparisons are 
appropriate when the fundamental equations are cast in the same interface, and that is possible at this moment only 
using Ashtekar's connection variables 18|. More precisely, the authors' recent studies showed the following: (a) 
the original set of dynamical equations proposed by Ashtekar already forms a weakly hyperbolic system |]l9|, (b) 
by requiring additional gauge conditions or adding constraints to the dynamical equations, we can obtain a strongly 
hyperbolic system , (c) by requiring additional gauge conditions and adding constraints to the dynamical equations, 
we can obtain a symmetric hyperbolic system |l9|j20[| , and finally (d) based on the above symmetric hyperbolic system, 
we can construct a set of dynamical systems which is robust against perturbative errors for constraints and reality 
conditions |2|] {aka. A-system [p2[). 

Based on the above results (a)-(c), we developed a numerical code which handles gravitational wave propagation 
in the plane symmetric spacetime. We performed the time evolutions using the above three levels of Ashtekar's 
dynamical equations together with the standard ADM equation. We compare these for accuracy and stability by 
monitoring the violation of the constraints. We also show the demonstrations of our A-system (above (d)) in the 
succeeding paper (Paper II) , together with new proposal for controlling the stability. 

It is worth remarking that this study is the first one which shows full numerical simulations of Lorcntzian spacetime 
using Ashtekar's connection variables. This research direction was suggested |2^ soon after Ashtekar completed his 
formulation, but has not yet been completed. Historically, an application to numerical relativity of the connection 
formulation was also suggested ||l8|,^ using Capovilla-Dell- Jacobson's version of the connection variables |2^ , which 
produce an direct relation to Newman-Penrose's ^s. However here we apply Ashtekar's original formulation, because 
we know how to treat its reality conditions in detail [ p6p7| ], and how they form hyperbolicities. We will also describe 
the basic numerical procedures in this paper. 

The outline of this paper is as follows: In the next section, we review the mathematical background of the hyperbolic 
formulation briefly and present our fundamental dynamical equations. In §111, we describe our numerical procedures. 
Our experiments are presented in and §0 and we summarize them in §VI. The Appendix |a| is for showing 
Ashtekar's basic equations in our notation, and we also present our experiments based on the Maxwell equation in 
the Appendix We also introduce briefly the discussion in our Paper II in the Appendix 



II. HYPERBOLIC FORMULATIONS 
A. definitions, properties, mathematical backgrounds 

We say that the system is a first-order ( quasi-linear) partial differential equation system^ if a certain set of (complex- 
valued) variables (a = 1, • • • , rt) forms 

dtUo,^M'Po.{u)diup+No.{u), (2.1) 

where M. (the characteristic matrix) and M are functions of u but do not include any derivatives of u. If the 
characteristic matrix is a Hermitian matrix, then we say ( |2.l| ) is a symmetric hyperbolic system. 

Writing the system in a hyperbolic form is the essential step in proving the system is well-posed. Here, well- 
posedness of the system means (1°) existence (of at least one solution u), (2°) uniqueness (i.e., at most solutions), 
and (3°) stability (or continuous dependence of solutions {u} on the Cauchy data). The Cauchy problem under 
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weak hyperbolicity is not, in general, C°° well-posed. The symmetric hyperbolic system gives us the energy integral 
inequalities which are the primary tools for studying the stability of the system. Well-posedness of the symmetric 
hyperbolic is guaranteed if the characteristic matrix is independent of u, while if it depends on u we have only the 
limited proofs for the well-posedness. From the mathematical point of view, proving well-posedness with less strict 
conditions is an old but active research problem. 

We can define another hyperbolic system between the weakly and symmetric levels. For example, we say we have 
a strongly hyperbolic (or diagonalizable hyperbolic system, if the characteristic matrix is diagonalizable and has 
all real eigenvalues. The inclusion relation is, then, 

symmetric hyperbolic G strongly hyperbolic G weakly hyperbolic, (2-2) 

(which means the symmetric hyperbolicity requires stronger conditions to be satisfied than the others). We do not 
repeat each level's features here (see §2 of [|l9|). However, at the strongly hyperbolic level, we can prove the finiteness 
of the energy norm if the characteristic matrix is independent of u (cf [^6|), that is one step definitely advanced than 
a weakly hyperbolic form. 

From the point of numerical applications, to write down the fundamental equation in an explicitly hyperbolic form 
is quite attractive, not only for its mathematically well-posed features. It is well known that a certain flux conservative 
hyperbolic systems of equations is taken as an essential formulation in the computational Newtonian hydrodynamics 
Ipsf . There is also an effort for implementing the boundary condition by using the characteristic speed (eigenvalues) 
of the system . 



B. Hyperbolic formulations of the Einstein equation 

As was discussed by Geroch | |30| , most physical systems can be expressed as symmetric hyperbolic systems. However, 
the standard ADM system does not form a first order hyperbolic system. This can be seen immediately from the fact 
that the ADM dynamical equations, 

da^j = -2NK,j + VjN, + V,N, , (2.3) 

+ (VjiV")i^„,, + (V,7V™)X,„, + N-^VmK,, , (2.4) 

have Ricci curvature ^^^Rij which involves second derivatives of the three-metric 7^ by definition. (The notation here 
is the standard one. Kij is the extrinsic curvature, N and iV' are the lapse and shift vector, respectively. V denotes 
a covariant derivative on the three-surface.) For our later convenience, we also write down the ADM constraint 
equations, 

(.ADM _ (3)^ ^ (^^^)2 „ J^^^J^^J _ ^2.5) 

^ADM^ ._ Vj (xy - ^HrK) « 0, (2.6) 

which are called the Hamiltonian and momentum constraint equations, respectively. 

So far, several first order hyperbolic systems of the Einstein equation have been proposed; some of them are 
symmetrizable (strongly hyperbolic) or symmetric hyperbolic systems ^2 14 3l|j3^ . There are many variations 



in the methods for constructing higher hyperbolic systems, but the number of fundamental dynamical variables is 
always results in larger than that of ADM (see a brief example by Anderson- York(99) in (ij]). Several numerical 
tests are reported (as we referred in the Introduction) using a particular hyperbolic formulation, but no numerical 
comparisons between these formulations are reported [Q. 

Using Ashtekar's formulation, we can compare three levels of hyperbolicity in the same interface (same fundamental 
variables) as we describe next. 
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C. Hyperbolic formulations in the Ashtekar formulations 



We present here our fundamental dynamical equations. Our notations and a more detailed review are presented in 
the Appendix but we repeat them here if necessary. 

The new basic variables are the densitized inverse triad, E^, and the S0(3,C) self-dual connection, Af, where the 
indices indicate the 3-spacetime, and a, 6, • • • are for S0(3) space. The total four-dimensional spacetime is 

described together with the gauge variables N,N^,Aq, which we call the densitized lapse function, shift vector and 
the triad lapse function. The system has three constraint equations, 

C^^'':^{t/2)e''\ElElF.^^0, (2.7) 
C^f := -Ft.Ei « 0, (2.8) 
= ^^^K - 0, (2.9) 



lASH 
--Ga 



which are called the Hamiltonian, momentum, and Gauss constraints equation, respectively. The dynamical equations 
for a set of {El,Af) are 

dtEi = -iV,{e^\ NEiEl) + 2V,{N^^Ej) + zA^eab' E^ (2.10) 
dtAI = -le'^K NElF^j + N^Ff^ + V,A^, (2.11) 

where F° := 2d[iA^^ — ie°'hcA\A'^j is the curvature 2-form. 

We have to consider the reality conditions when we use this formalism to describe the classical Lorentzian spacetime. 



As we review in § |A 2| , the metric will remain on its real- valued constraint surface during time evolution automatically 



if we prepare initial data which satisfies the reality condition. More practically, we further require that triad is real- 



valued. But again this reality condition appears as a gauge restriction on ^Iq, (All), which can be imposed at every 
time step. In our actual simulation, we prepare our initial data using the standard ADM approach, so that we have 
no difficulties in maintaining these reality conditions. 



The set of dynamical equations (2.10) and (2.11) [hereafter we call these the original equations] does have a weakly 
hyperbolic form Jl9| , so that we regard the mathematical structure of the original equations as one step advanced from 
the standard ADM. Further, we can construct higher levels of hyperbolic systems by restricting the gauge condition 
and/or by adding constraint terms, C^^^,C^j^^ and C^^^, to the original equations ||l^. We extract only the final 
expressions here. 



In order to obtain a symmetric hyperbolic system q we add constraint terms to the right-hand-side of (2.10) and 



(2.11). The adjusted dynamical equations, 

dtK = -I'Dji^^'a NEiEl) + 2V,{N^'Ej) + lA^tab^ El + P\bCf'^\ (2.12) 

where P\i, = N'Sab + iNeab^K, 

dtA"} = -ie'^ NElF^^^ + N^Pp^ + V,Al + Q,"C^^" + R^" C^//, (2.13) 

where ee e-^NE^, R^" = ie-^Ne^'bE^Ei 

form a symmetric hypcrbolicity if we further require the gauge conditions, 

A"^ = A'^N\ d,N = 0. (2.14) 

We remark that the adjusted coefficients, P'^ab,Qi,Ri-'°', for constructing the symmetric hyperbolic system are 
uniquely determined, and there are no other additional terms (say, no C^^^,C^/^ for 9fi?*, no Cq^^ for dtAf) 



9[. The gauge conditions, (2.14), are consequences of the consistency with (triad) reality conditions. 



Iriondo et al |34] presented a symmetric hyperbolic expression in a different form. The differences between ours and theirs 
are discussed in MJ20 
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We can also construct a strongly (or diagonalizable) hyperbolic system by restricting to a gauge N'- ^ 0, ±Ny^ 

Or we can 



(where 7" is the three-metric and we do not sum indices here) for the original equations (2.10), (2.11) 
also construct from the adjusted equations, ( ^.12 ) and ( 2.13 ), together with the gauge condition 



As for the strongly hyperbolic system, we hereafter take the latter expression. 

In Table |, we have summarized the equations to be used throughout the remainder of this article. 



(2.15) 



system 


variables 


Eqs of motion 


remark 


ADM 




( 


2.3), (2.4) 


"standard ADM" 


I Ashtekar (weakly hyp.) 


{El, A?) 


( 


2.1c 


), 


2.11 


) (original) 


"original Ashtekar" 


II Ashtekar (strongly hyp.) 




( 


2.1s 


), 


2.1; 


) (adjusted) 


( 


2.15 


) required 


Ill Ashtekar (symmetric hyp.) 


[El, At) 


( 


2.1s 


), 


2.1J 


) (adjusted) 


( 


2.14 


) required 


(in [V) Ashtekar (adjusted) 


[El, At) 


( 


5.1) 


, (" 


.2) (adjusted with n) 





TABLE I. List of systems that we compare in this article. 



III. NUMERICAL METHOD 
A. Overview 

We coded up the program so as to compare the evolutions of spacetime with different set of dynamical equations 
but with the common conditions: the same initial data, the same boundary conditions, the same slicing condition and 
the same evolution scheme. 

We consider the plane symmetric vacuum spacetime without cosmological constant. This spacetime has the true 
freedom of gravitational waves of two polarized (+ and x) modes. We apply the periodic boundary conditions to 
remove any difficulties caused by numerical treatment of the boundary conditions. The initial data are given by 
solving constraint equations in ADM variables, using the standard conformal approach by York and O'Murchadha 
p5[ . When we use Ashtekar's variables for evolution, we transform the ADM initial data in terms of Ashtekar's 
variables. The results are analyzed by monitoring the violation of the constraint equations which are expressed using 
the same (or transformed if necessary) variables. 

We describe our procedures in the following subsections in detail. 



B. metric and the initial data construction 



We consider the plane symmetric metric, 

ds^ = (-A^^ + NxN'-')dt^ + 2Nxdxdt + ^xxdx^ + ^yydy^ + 7^^dz^ + 2^y^dydz (3.1) 

where the components are the function of N{x,t), Nx{x,t),^xx{x,t),^yy{x,t),^zz{x,t),^yz{x,t). N and A^ are called 
the lapse function and the shift vector. 



We prepare our initial data by solving the ADM constraint equations, (2.5) and ( |2.6| ), using the conformal approach 
p5[ . Since we consider only the vacuum spacetime, the input quantities are the initial guess of the 3-metric 7^, the 
trace part of the extrinsic curvature trAT, and the transverse traceless part of the extrinsic curvature Att- For 
simplicity, we impose Att — and trA" = Kq (constant). The Hamiltonian constraint, then, becomes an equation 
for the conformal factor, -0: 

SAV- 8l=d,{f^ ^dj4,) - i?V + \{Ko)^i^\ (3.2) 
V7 3 
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where 7 = det ^ij. The momentum constraint is automatically satisfied by assumption. The initial dynamical 
quantities 7ij , Kij are given by the conformal transformation, 

7y = '0'*7jj, = ^^p'^%Ko. (3.3) 



We solve (3.2) under the periodic boundary conditions using the incomplete Cholesky conjugate gradient (ICCG) 
method. 

We should remark here that we have to assume non-zero Kq for a model of gravitational pulse waves under the 
periodic boundary conditions in this plane symmetric spacetime. This can be seen as follows. Suppose we set Kq = 0. 



From (3.2), we get 



= ^ J ^R^dx. (3.4) 



If we set the boundary as x = [A, B] and impose the periodic boundary conditions, then eq. (3.4) becomes 



dx'^\x=A - dx^\x=B ^ [ —j^-rz] I ^/^Ripdx. (3.5) 

VV75 Jx=a=bJb 

However when there exist a gravitational wave pulse which produces i? ^ in the region, this equation gives dx^\x=A = 
dx4'\x=B, which is inconsistent with the periodic boundary conditions. Therefore we need to assume non-zero Kq in 
order to compensate the curvature which is produced by the pulse waves. 

Actually the trace of the extrinsic curvature appears only in the quadratic form, so we can interpret that our 
(background) spacetime is either expanding, Kq < 0, or contracting, Kq > 0. However, this fact indicates that there 
is no known exact solution to compare with. If the background aspacetime is allowed to be flat (Kq = 0), then we know 
there is a series of exact solutions which describes a collision of plane gravitational waves which were originally found 
by Szekerez and Khan-Penrose [^6| . The formation of a curvature singularity after such colliding waves is known to be 
generic, but that is not generalized to the expanding background (as discussed using numerical simulations p7| , ^ ). 

We can set two different modes of gravitational waves. One is the -I— mode waves, which is given by setting a 
conformal guess metric as (in a matrix form) 

/I \ 

lij — sym. 1 -I- aexp(— 6(2:: — c)^) (3.6) 

\sym. sym. 1 — a cxp(— 6(a; — c)^) / 

where a, 6, c are parameters. The other is the x-mode waves, given by 

/ 1 \ 

7y = I sym. 1 aexp(-5(a; - c)^) 1 (3.7) 
\ sym. sym. 1 / 

where a, 6, c are parameters again. Both cases, we expect non-linear behavior when wave's curvature becomes quite 
large compared to the background. In the collision of a -I— mode wave and a x-mode wave, we also expect to see the 
mode- mixing phenomena which is known as gravitational Faraday effect. These effects are confirmed in our numerical 
simulations. 

C. Transformation of variables: From ADM to Ashtekar 

We need to transform the dynamical variables on the initial data when we evolve them in the connection variables. 
We list the procedure to obtain (i?^, ^°) from {-fij,Kij). This procedure is used also when we evaluate the constraints, 
Cff^^'^Mi^'^cf" foi' the data evolved using ADM variables. 

From the three- metric 7^^ to i?* : 
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1 . Define the triad corresponding to tlie three-metric 7^ . We take 



El El 

r?2 rp2 r?2 

Ex E tj^ 

El eI El 



/l^ 

622 623 
632 633 



and set simply 623 = 632. The relation between the metric and the triad becomes 



2 I 2 _ 

^22 + 633 — ^yy. 



-23 



+ 633 = 7^27 (622 + 633)623=71^2 



(3.8) 



(3.9) 



For the case of +-mode waves, we define natm'ally, 622 — \Jlyy^ 633 — ^/Yzz, 623 — 0. For x-mode waves, we also 
take a natural set of definitions, 622 = 633 = [(7^^ + (7^^, - 7^2)^/^)/2]^/^ and 623 = 7yz/2622 which are given 



by solving 632 + 633 = ^yy and 2622623 = lyz- 

2. obtain the inverse triad from triad Ef. 

3. calculate the density, 6, as 6 = detii^f. 

4. obtain the densitized triad, E'* — eE^. 
From three- metric {jij,Kij) to Af: 

1. prepare the triad Ef and its inverse E^. 

2. calculate the connection 1-form uj'^'^ — E^'-'-^iE'^j^. This is expressed only using partial derivatives as | 



E^'di^E^^ - E^dE'^'E^^d^kE^^ 



(3.10) 



3. At = -K.jW" - |e%c^f. 



D. Transformation of variables: From Ashtelcar to ADM 



In contrast to the previous transformation, we also need to obtain (jij^Kij) from {El,A^) when we evaluate the 
metric output or ADM constraints when we evolve the spacetime using connection variables. This process is only 
required at an evaluation times, not required at every time step (unless we use the gauge condition which is primarily 
defined using ADM quantities). 

From densitized inverse triad E^^ to three-metric 7^ : 

1. calculate the density e as e = {det E^Y^^ . 

2. get the three inverse metric as 7'^ = El^Ei/e^ . 

3. obtain 7^ . 

From (i?^,X°) to the extrinsic curvature Kij-. 



E^''ViE^ and cj"'"' := E^^co"", and a relation 

Q [a6c] r\ [bc]a a [be] . 6[ca] . c[a6] abc . cba abc 



Using the densitized triad, eq. ( [3.1C| ) can be also expressed as 

= ^E^'i^i^E]]) + ^&'mK{d,El) + -^E,a&'Ei{djEl), taking [he]. 
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1. prepare the un-densitized inverse triad, i?* = E'a/^- 

2. prepare triad Ef . 

3. calculate the connection l-form e°- ^c'-Jl'^ ■ 

4. calculate Zf, which is defined as Zf -ylf + ^e'^^uo^^i^ K.jE^"), and get = i^ja- 

E. Gauge conditions 

We evolve the initial data with different evolution equations and compare its accuracy/stability. As we summarized 
in Table I, we will compare time evolutions between ADM and Ashtekar (of the original system I) in § [V A , and three 



of Ashtekar's systems (I, II and III: weakly, stronly and symmetric) in § IVB . We, then, consider alternative system 
(adjusted k system) in §0 

Here we comment again on our choice of the slicing (gauge) condition. As for the primary tests of this subject, we 
apply the simplest slicing conditions we can take. That is, 

(1) the simplest geodesic slicing condition for the lapse function, 

(2) the simplest zero shift vector A^^ = 0, and 



(3) the natural choice of triad lapse function Af, = AfN^ [= if A''^ = 0, which is suggested from (2.14) or (2.15)] 



However, in the Ashtekar formalism, the densitized lapse function A^ is the fundamental gauge quantity (rather than 
N). Therefore we try two conditions for the lapse, 

(la) the standard geodesic slicing condition A^ = 1, which will be transformed to A^ = 1/e when we apply this 
condition in Ashtekar's evolution system, and 

(lb) the densitized geodesic slicing condition A^ — 1, which will be transformed to A^ = e when we evolve the system 
using ADM equations. 

In practice, such a transformation using the density e will not guarantee that the Courant condition holds if we fix 
the time evolution step At ^ Therefore we need to rescale the transformed lapse [ A^ in (la), A^ in (lb)] so that it 
has a maximum value of unity, in order to keep our evolution system stable. 

If we apply the standard geodesic slice, then we can compare the weakly hyperbolic system with the symmetric hy- 
perbolic one. Similarly if we apply the densitized geodesic slice, then we can compare the (original) weakly hyperbolic 
system with the strongly hyperbolic one. 

F. Time integrating scheme 

We applied two second order evolution schemes, and confirmed that they give us nearly identical results. 
One is the so-called iterative Crank-Nicholson scheme (cf. ||39[] ), which is now becoming the standard in the numerical 
relativity community. Suppose we have a dynamical equation in the form 

dtu{x, t) = f{u{x, t), dxu{x, t)). (3.11) 



■* This is from tfie original definition of A?, At — t<J°" - (i/2)e"bc t^^. 

^We here remind the reader of the stability condition, N dd < Aa; for a standard forward-time centered-space (FTCS) scheme 
for a simple wave equation, in a (At, Aa;)-spaced numerical grid. Note that this condition will be changed due to the choice of 
the evolution scheme and the equations of the system. 
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Then the scheme for updating m at a point x from t to t + At consists of the following steps. (1) use data ont = t for 
right-hand-side and update u{x, t) for At step as u{x, t + At), 

u{x,t + At)-u{x,t) 

^ f{u{x,t),d^u{x,t)). (3.12) 

(2) take the average of u{x, t) and u{x, t + At), and let it represent a half-step value, (say v{x, t + At/2)). (3) update 
u{x, t) for At step again using u{x, t + At/2) in the argument of the right-hand-side, 

''i^^t + ^^^-^i^'i) ^ t At/2),dMx, t + At/2)). (3.13) 

(4) perform the above (2) and (3) steps once more (we suppose two- iteration Crank-Nicholson scheme), and take 
u{x, t + At) to be the evolved quantity. 

The other scheme we applied is the Brailovskaya integration scheme, which is a second order predictor-corrector 



method Q| and rather easy to code. The first step (predictor step) is the same as (3.12), and the second step 
(corrector step) simply switches the right-hand-side using the updated u(x, t + At) to be 

u(a;, t -t- At) — u(a;, t) , . x ~/ . xx /„ ,x 
= /(u(x, t + At), dM^, t -f At)). (3.14) 

Note that all derivatives here in the right-hand-side are assumed to use central difference. 

The latter scheme is quite simple, but gives us reasonably accurate and stable evolutions for our problems. We 
confirmed that both give us nearly identical evolutions (which will be shown in Fig|| (b)), but the Brailovskaya 
method requires less computational time. 

G. Checking the constraints 

We compare the violation of the constraint equations during the time evolution. We have ADM constraint equations, 



C^DM g^^^ qAdm j^j^j (m)]^ g^^^ j^lgQ Ashtekar's constraint equations, C^^h^^ash qASH [(|||)^ Q ^^^^ 



(2.9), respectively] . By means of the transformation between (7^ , Kij ) and {E^ , A"; ) , we can evaluate ADM constraints 
even if we evolved the system using Ashtekar's variables and vice versa. 

We measure the violation of a constraint by its (i) maximum, max2;|C(a;)|, (ii) LI norm, (X]"=i ^^'^ (iii) 

L2 norm, (X]"=i |C(a;)p/na;)"'^^^, where is the number of grid points. When we compare them during the evolution, 
we measure them at the same proper time, r, for the two different evolution systems. The proper time is defined 
locally as dr = Ndt, which is also dr = eNdt, but here we apply its averaged value on the whole t =constant surface, 

(say(iV)^(lK)E:=iA^), 

T = I {N)dt, or T= / {eN)dt, (3.15) 
Jo Jo 

to characterize the "time" of evolution. 

The numerical code passed convergence tests, and the results shown in this article are all obtained with acceptable 
accuracy. In Fig.|l|, we show a result of convergence tests. We show the convergence behaviour of our initial data solver 
in Fig.0(a), together with the convergence behaviour of the evolution codes both for ADM and Ashtekar variables in 
Fig.|l|(b) and (c). We plotted the residual of the Hamiltonian constraint solver when it was minimized, and the L2 
norm of C^^^ and C^^^ for the evolution of +-mode single pulse wave. (The model is described in the next section.) 
We can see all errors are diminished by finer resolutions. The order of the convergenc^ is 1.98 to 2.01 for the initial 
data solver, and at best 1.96 (e.g. at r = 0.5) for the evolution code. 



Here we used the definition of the order of convergence following Bona et. al. 
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All the results we present in this article are obtained using 401 grid points for the range x = [—5, +5]; that is 
gravitational waves traverse the entire numerical region in proper time 10 if the background expansion Kq is close to 
zero. We use the Courant number v = At/ Ax — 0.2. 

We coded all our fundamental quantities (metric, gauge variables, • • •) as complex, but we observed that the 



evolution from our initial data never violate its metric reality conditions. Due to the gauge condition for Aq, (All), 
we also confirmed that our evolutions preserve the triad reality condition. 

IV. EXPERIMENTS 1: DIFFERENCES BETWEEN HYPERBOLICITIES 

In this section, we examine the accuracy/stability of the numerical evolutions comparing the different hyperbolic 
systems. We begin showing how the evolution of Ashtekar's equations look, comparing with those of the ADM 
equations. 

A. ADM vs Ashtekar 

We start by describing our model, plane wave propagation in expanding/collapsing spacetime. We prepare the 
initial data with one or two gravitational pulse waves in our numerical region. The pulses then start propagating in 
both ±x directions with the light speed, and appear on the other side of the numerical region due to the periodic 
boundary condition. When the pulses collide, then the amplitude seems simply to double, as they are superposed, 
and the pulses keep traveling in their original propagation direction. That is, we observe something like solitonic wave 
pulse propagation. 



As we mentioned in §111 B, we have to assume our background not to be flat, therefore there is no exact solutions. 
The reader might think that if we set |tr if | to be small and pulse wave shapes to be quite sharp then our simulations 
will be close to the analytic colliding plane wave solutions which produce the curvature singularity. However, from 
the numerical side, these two requirements are contradictory (e.g. sharp wave input produces large curvature which 
should be compensated by |tr K\ in order to construct our initial data). Thus it is not so surprising that our waves 
propagate like solitons, not forming a singularity. 

In Fig.|| (a), we plot an image of wave propagation (metric component gyy) up to t = 10, of +-mode pulse waves 
initially located at x = ±2.5. We took a small negative Kq, so that the background spacetime is slowly expanding. 

Fig.^ (b), then, tells us that our ADM evolution code and Ashtekar's variable code give us identical evolutions. We 
plotted a snapshot of gyy on the initial data (which is common to all models here), and its snapshot at r = 10.0. The 
fact that all four lines (ADM/ Ashtekar, of their Brailovskaya/Crank-Nicholson evolution schemes) overlapped clearly 
indicates that we are showing exact evolutions. 

We also plotted a typical evolution of the fundamental dynamical quantities and Ay in Figs.|| (c) and (d). 

We next compare constraint violations by Ashtekar's equation with that of ADM. In FigJ^, we plot the L2 norm of 
C^^^ and C^^^. We see that ADM evolution shows less violation in measuring C^^^, and the Ashtekar evolution 
shows less violation in measuring CASH. The magnitudes of these violations are similar. Thus, we believe, these 
violations are within the numerical truncation errors in the process of numerical transformation of variables (ADM 
to Ashtekar/ Ashtekar to ADM), and therefore it is not appropriate to conclude here which formulation is better. 

As the reader may guess, the violations of constraints reduce if the background spacetime is expanding {Kq < 0). 
Therefore we will use the collapsing background spacetime {Kq > 0) hereafter for presentations, with the expectation 
of having more non-linear effects; however this direction also stops the evolution after the finite time. (e.g. for the 
flat initial data with Kq = +0.025, the spacetime will collapse to zero volume around t — 60.) 

B. Comparison between hyperbolicities 

We here present our comparisons of the accuracy and/or stability between the different hyperbolicities. Since all 
examples we show in this section are not the case of unstable evolution (no exponential growth of constraint violation) , 
our experiments can be said as the comparisons of the accuracy of the evolution, conservatively. 
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We first compare the (original) weakly hyperbolic system [system I in Table || with the strongly hyperbolic system 
[system II in Table ^. This comparison can be done under the densitized geodesic slicing condition, N = I. We 

prepare two initial gravitational pulses (both +-+ or x-x modes) and take our background spacetime to be collapsing 
{Kq > 0). In Fig.^, we show the constraint errors, C^^^ and C^/^. In both two situations, we observe that the 
strongly hyperbolic system has slightly improved the violation of the constraints, but we can not see the orders of 
magnitude differences. 

Similarly, we next compare the (original) weakly hyperbolic system [system I in Table |] with the symmetric 
hyperbolic system [system III in Table ||] . This comparison can be done under the standard geodesic slicing condition, 
A'' = 1. We repeat the same experiments as above and show plots in Fig.^. We again see that the symmetric hyperbolic 
system slightly improves the situation, but not so drastically. 

From both Figs. ^ and ^ we see that the strongly and symmetric hyperbolic systems produce less violation of 
constraints than the original weakly hyperbolic system. Therefore one conclusion is that adjusting the equation of 
motion with constraint terms does definitely make the system accurate. However the constraint violation remains the 
same order of magnitude. 

From each figures, we may conclude that higher level hyperbolic system gives us slightly accurate evolutions. 
However, if we evaluate the magnitude of L2 norms, then we also conclude that there is no measurable differences 
between strongly and symmetric hyperbolicities. This last fact will be supported more affirmatively in the next 
experiment. 



V. EXPERIMENTS 2: ANOTHER WAY TO CONTROL THE ACCURACY/STABILITY 

The results we have presented in the previous section indicate that both strongly and symmetric hyperbolic systems 
show better performance than the original weakly hyperbolic system. These systems are obtained by adding constraint 



terms (or "adjusted" terms) to the right-hand-side of the original equations, ( 2.10 ) and ( 2.11 ). In this section, we 
report on simple experiments in changing the magnitude of the multipliers of such adjusted terms. 

We consider the following system, where the equations of motion are adjusted in the same way as before, but with 
a real-valued constant multiplier k: 

dtK = -i'D.ie^K. NEiEl) + 2V,{N^^Ej) -f iA^eab' K + KP\tC^^''', (5.1) 

where P\b = N'Sab + iNeab^El, 

dtA"} ^ -ie^-^ NE'iF^j + N'F^^ + V.A^^ + nQ^^Cf^ + Ki?.^" C^/f , (5.2) 

where Q'^ = e-^NEf, R,^" = ie~^ Ne^^bEtsl. 



The set of ( f3.l| ) and (5.2) becomes the original weakly hyperbolic system if k 0, becomes the symmetric hyperbolic 
system if k = 1 and N = const., and remains strongly hyperbolic systems for other choices of k except k = 1/2 
which only forms weakly hyperbolic system. We remark again that the coefficients for constructing the symmetric 
hyperbolic system are uniquely determined. 

We tried the same evolutions as in the previous section for different value of n. In Fig.^, we plot the L2 norm of 
the Hamiltonian and momentum constraint equations, and C^f^- We checked first that k — Q and 1 produce 

the same results as those of weakly and symmetric hyperbolic systems. What is interesting is the case of k = 2 
and 3. These ks produce better performance than the symmetric hyperbolic system, although these cases are of 
strongly hyperbolic levels. Therefore, as far as monitoring the violation of the constraints is concerned, we may say 
the symmetric hyperbolic form is not always the best. We remark that the negative k will produce unstable evolution 
as we plotted, while too large positive k will also results in unstable evolution in the end (see k = 3 lines). 

We also tried similar experiments with the vacuum Maxwell equation. The original Maxwell equation has symmetric 
hyperbolicity, and additional constraint terms (with multiplier k) reduce hyperbolicity to the strong or weak level. 
We show the details and a figure in the Appendix but in short there may be no measurable differences between 
strongly and symmetric hyperbolicities. 
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These experiments in changing k are now reported in our Paper II |^ more extensively. There, we propose 
a plausible explanation why such adjusted terms work for stabilizing the system. We introduce the idea in the 
Appendix |^. Briefly, we will conjecture a criterion using the eigenvalues of 'adjusted version' of the constraint 
propagation equations. This analysis may explain the appearance of phase differences between two systems, which is 
observed in Figs.^, || and ||. 



VI. DISCUSSION 



Motivated by many recent proposals for hyperbolic formulations of the Einstein equation, we studied numerically 
these accuracy/stability properties with the purpose of comparing three mathematical levels of hyperbolicity: weakly 
hyperbolic, strongly hyperbolic, and symmetric hyperbolic systems. We apply Ashtekar's connection formulation, 
because this is the only known system in which we can compare three hyperbolic levels with the same interface. 

Our numerical code demonstrates gravitational wave propagation in plane symmetric spacetime, and we compare 
the "accuracy" and/or "stability" by monitoring the violation of the constraints. Actually, our experiments in §rV 
were the comparisons of accuracy in evolutions, while in §^ we observed cases of unstable evolutions. By comparing 
with the results obtained from the weakly hyperbolic system, we observe the strongly and symmetric hyperbolic 
system show better properties with little differences between them. Therefore we may conclude that higher levels of 
hyperbolic formulations help the numerics more, though its differences are small. 

However, we also found that the symmetric hyperbolic system is not always the best for controlling accuracy or 
stability, by introducing a multiplier for adjusted terms in the equations of motion. This result suggests that a certain 
kind of hyperbolicity is enough to control the violation of constraint equation. In our case it is the strongly hyperbolic 
level. This statement is supported by an experiment in Maxwell system as we describe in Appendix 

The remaining question is: why we can get the better performance by adding constraint terms in the dynamical 
equations? The added terms are basically error terms during the evolution for its original dynamical equations. 
Nevertheless, these terms improve the accuracy of the evolution. We now have a plausible way to explain the reason 
which is discussed in our Paper II ji^ (its brief introduction is in Appendix ^ in this article). There we evaluate 
the eigenvalues of the adjusted version of the constraint propagation equations, and propose a criteria for obtaining 
stability of the system. In some cases, for example, the decay/growth of the constraints can be predicted the signature 
of the eigenvalues of the adjusted version of the constraint propagation equations. In |43|, we will discuss this point 
in detail together with a numerical demonstration of A-systems [^,^. There we also show that some choices of 
adjusted terms may produce unstable evolution. 

To conclude, we are glad to announce that Ashtekar's connection variables have finally been applied in numerical 
simulations. This new approach, we hope, will contribute to understanding further of gravitational physics, and will 
open a new window for peeling off interesting non-linear natures together with a step to numerical treatment of 
quantum gravity. 
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APPENDIX A: ASHTEKAR'S FORMULATION OF GENERAL RELATIVITY 



We give a brief review of the Ashtekar formulation and the way of handling reality conditions. This appendix is for 
describing our notations. 
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1. Variables and Equations 



The key feature of Ashtekar's formulation of general relativity |17| is the introduction of a self-dual connection as 
one of the basic dynamical variables. Let us write the metric using the tetrad E'^ as ~ E^^E^r\ij Define 
its inverse, iJj, by E*^ :— E^g^^'^rju and we impose = as the gauge condition. We define S0(3,C) self-dual 
and anti self-dual connections '^A't '■= T (*/2)e"bc wj^'^, where uj^J is a spin connection 1-form (Ricci connection), 



:— E^"\I uE'l- Ashtekar's plan is to use only the self-dual part of the connection and to use its spatial part 



as a dynamical variable. Hereafter, we simply denote as ^J^. 

The lapse function, N , and shift vector, A^', both of which we treat as real- valued functions, are expressed as i?Q = 
(1/A'^, —N'^/N). This allows us to think of Eq as a normal vector field to S spanned by the condition t = =const., 
which plays the same role as that of Arnowitt-Deser-Misner (ADM) formulation. Ashtekar treated the set (E'*, A'^) 
as basic dynamical variables, where E^^ is an inverse of the densitized triad defined by i?* := ei?*, where e := detEf 
is a density. This pair forms the canonical set. 

In the case of pure gravitational spacetime, the Hilbert action takes the form 

S = j d''x[{dtA'i)K + i^imKElF^,^''''c - e^ATV - N^F^^Ei + A^ V^E^], (Al) 

where N := e^^N, F^^ :— 2d[^A'l^~ie°'bc A^^A^ is the curvature 2-form, A is the cosmological constant, ViEl := diE^— 
ieah''A\Ei, and = deti^^ = (det£;f)2 is defined to be det£;^ = {l/6)€'''"'ujkKK^c> where e^fe := CabcE^E^E^ 
and ei-jk ■= e'^e^fc |. 

Varying the action with respect to the non-dynamical variables N, and Aq yields the constraint equations, 

C^SH (i/2)e'^'', ^a^fci^j - A det ^ w 0, (A2) 
CtF^ ~Fr,Ei « 0, (A3) 
C^f T^^El « 0. (A4) 

The equations of motion for the dynamical variables (i?^ and A°;) are 

dtK = -*2?,(e"'a NEiEl) + 2V,{N^^E^) + lA^'^eab" K, (A5) 

dtAf = -^e'^^ NElF^''^ + F^^ + V.A^ + KNEf , (A6) 

where P^Xf := d^Xl' - ieahA)Xl\ for X]^ + X^' = 0. 



2. Reality conditions 

In order to construct the metric from the variables {El^, A"; , N, N"^), we first prepare the tetrad Ej as Eq = 

{l/eN, -N'/eN) and E^ = (0,^^/e). Using them, we obtain the metric gf"" such that gf"" := E'^E^'j-q". 

This metric, in general, is not real-valued in the Ashtekar formulation. To ensure that the metric is real-valued, we 
need to impose real lapse and shift vectors together with two metric reality conditions; 

Im(^:^^'^) = 0, (A7) 

W'^ Reie'''"'E^El"VkEi'^) = 0, (A8) 



'^We use — 0, ■ ■ ■ ,3 and i, j = 1, • • ■ , 3 as spacetime indices, while I,J = (0), ■ ■ ■ , (3) and a,b = (1), ■ ■ ■ , (3) are 50(1, 3), 
SO(3) indices respectively. We raise and lower - ■ ■ hy g^" and gpi^ (the Lorentzian metric); /, J, • ■ ■ by 77" = diag(-l, 1,1,1) 
and r)ij\ • ■ ■ by j'^ and ^ytj (the three- metric); a, 6, ■ ■ ■ by 5"'' and Sab- We also use volume forms eabc- eabce'''"^ = 3!. 

*When {i,j,k) = (1,2,3), we have e^^fc = e, e.jk = 1, e'^'' = e~'^, and l'^'' = 1. 
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where the latter comes from the secondary condition of reality of the metric Im{9t(i?* = and we assume 

deti^ > (see 

For later convenience, we also prepare stronger reahty conditions, triad reality conditions. The primary and sec- 
ondary conditions are written respectively as 



^ n 



and 



:= lm(^^) = 0, 
liaidtEa) = 0. 



(A9) 
(AlO) 



Using the equations of motion of E^, the gauge constraint 
primary condition 



the metric reality conditions ( |A7D , ( [A8| ) and the 

(All) 



we see that (AlO) is equivalent to 

Re{A^) = d^{N)E'"- + {ll2e)ElNE^''djEl + iV'Re(^n, 
or with un-dcnsitized variables, 

Re(ytg) = d,{N)E"' + N' Yie{A1). (A12) 

From this expression we see that the secondary triad reality condition restricts the three components of the "triad 
lapse" vector ^q. Therefore (All) is not a restriction on the dynamical variables (i?* and A";) but on the slicing, 
which we should impose on each hypersurface. 

Throughout the discussion in this article, we assume that the initial data of {El^^AD for evolution are solved so as 
to satisfy all three constraint equations and the metric reality condition (A7) and (|A8|). Practically, this is obtained, 
for example, by solving ADM constraints and by transforming a set of initial data to Ashtekar's notation. 



APPENDIX B: EXPERIMENTS USING THE MAXWELL EQUATION 



In this appendix, using the Maxwell equation of the vacuum field, we show that the symmetric hyperbolic system 
does not change the stability feature drastically. The result here supports the discussion in More detail analysis 
can be found in our Paper 11 . 

The Maxwell equation has two constraint equations. 



Ce := d,E' « 0, Cb := d,B' 



0, 



and two dynamical equations 



for the field {Ei.Bi). 



Suppose we have adjusted (B2) using the constraint terms, (Bl), with a multiplier, k. 

dtE, = ce^^'djBk + k,Ce. dtB, - -ce^^^'djEu + k,Cb 
where = (k, k, k) for simplicity. This matrix expression 



dt 



-ce 



di 



(Bl) 



(B2) 



(B3) 



(B4) 



immediately tells us its hyperbolicity depending on k as follows: The system, (B4), becomes symmetric hyperbolic 
form when k = (that is the original Maxwell equation), becomes weakly hyperbolic form when k = ±c, and becomes 
strongly hyperbolic otherwise. The eigenvalues of the dynamical equation can be written as (c, c, — c, — c, k, k). 
We made a numerical code to demonstrate a propagation of plane electro-magnetic wave. 



i;Xx,t) = (0,0,-i=sin(^-ct)), 

B\x,t) = (-isin(^ - ct), isin(^^ - ci),0) 



(B5) 
(B6) 
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in 2-dimensional spacetime with periodic boundary condition. We use ( |B(]| ) as our initial data, and monitor its 
numerical error during its evolution by evaluating constraint equations and by checking error from the exact solution. 
The error itself is quite small, but as we show in Fig.^ we found the difference due to the multiplier of the adjusted 
terms k. We see that the symmetric hyperbolic equation show the best performance for the stability, but not showing 
so much different performance from the strongly hyperbolic system. 



APPENDIX C: WHY ADJUSTED EQUATIONS HAVE BETTER PERFORMANCE? 



Here, we try to explain briefly why the adjusted equations [( |5.1| ) and (5.2) for Ashtekar's system, (B3) for Maxwell 
equations] reduce the violation of constraints in the evolution. The detail explanations and numerical experiments 
are in our paper 11 |43| , and this appendix describes the essential idea of the mechanism. 

Suppose we have constraint equations, Ci « 0, C2 ~ 0, • • •, in a system. We, normally, monitor the error of the 
evolution by evaluating these constraint equations on the each constant-time hypersurface. Such monitoring, on the 
other hand, can be performed also by checking the evolution equations of the constraint, which we denote constraint 
propagation equations (cf. |^^ ). We, therefore, consider constraint propagation equation of which we transformed in 
Fourier components, C, 



(CI) 



The idea here is to estimate the eigenvalues of the matrix, M, after we took its leading order quantitiy in linealization 
against a particular background. Clearly, if all the real part of the eigenvalues are negative, then all constraints decays 
to zero along to the system's evolution. In our paper II we show that such a case can be obtained by adding 
'adjusted terms' both for Ashtekar's and Maxwell's systems. There we also show examples of unstable evolution by 
choosing adjusted terms which produce positive eigenvalues of M . The imaginary part of the eigenvalues are also 
supposed to contribute the appearance of the phase differences of the system. 

In this point, we can say that adjusted terms are responsible for obtaining the stable and/or accurate evolution 
system, and this is a way to control the stability of simulation, which effects more than the system's hyperbolicity. 
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FIG. 1. Examples of the convergence tests, (a) Convergence of the initial data solver [Hamiltonian constraint (3.2) solver]. 
We plot the residual of the conformal factor, L2 norm of \{^n — ipn-i) /'4>n\^ , when it converges, wh ere n is the iteration number 
in ICCG routine. The horizontal value is the amp litude of the gravitational wave [a in eq.(3.6)] where we assume +-mode 
single pulse wave and fix fe = 2.0, c = 0.0 in eq.( |3.6| ), and Ko = —0.025. According to the resolutions (grid points = 101, ■ ■ ■, 
801 for the range of a; = [—5, +5]), we see second order convergence, (b) Convergence behavior of ADM evolution code. L2 
norms of is plotted for the above initial data for the amplitude a = 0.2. We applied geodesic slicing condition (A'' = 1). 

(c) Convergence behavior of Ashtekar evolution code. L2 norms of Cf^^ is plotted for the above initial data for the amplitude 
a = 0.2. We applied geodesic slicing condition {N = 1). We can see clearly that the error norms in evolutions will decrease in 
high resolution cases. 
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FIG. 2. Images of gravitational wave propagation and comparisons of dynamical behaviour of Ashtekar's vari ables and 
ADM variables. We applied the same initial data of two +-mode pulse waves (a = 0.2,6 = 2.0, c — ±2.5 in eq.(3.6) and 
Ko — —0.025), and the same slicing condition, the standard geodesic slicing condition (TV = 1). Fig. (a) [top left] is the image 
of three-metric component gyy of a function of proper time r and coordinate x. This behaviour can be seen identically both in 
ADM and Ashtekar evolutions, and both with Brailovskaya and Crank-Nicholson time integrating scheme. Fig.(b) explains this 
fact by comparing the snapshot of gyy at the same propertime slice (r — 10), where four lines at r = 10 are looked identically. 
Figs.(c) and (d) [botoms] are of the real part of the densitized triad and the real part of the connection Ay, respectively, 
obtained from Ashtekar variables' evolution. 



18 




FIG. 3. Comparisons of the constraint violation by Ashtekar's equation with that of ADM. (a) L2 norm of the Ashtekar's 
Hamiltonian constraint equation, Cj^^ as a function of averaged proper time, (b) L2 norm of the ADM's Hamiltonian constraint 
equation, C^^^ as a function of averaged proper time. The plots are of the same parameters with those of Fig.0. 
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plus-mode wave propagation 




plus-mode wave propagation 




(bl) 
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cross-mode wave propagation 




(b2) 
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FIG. 4. Comparisons of the strongly hyperbolic system (Ashtekar II) with the weakly hyperbolic system (Ashtekar original) 
(TV = 1 slice). Figs. (a)s are of +-mode waves {a = 0.2,6 = 2.0, c = ±2.5 in eq.(3.G), while (b)s are of x-mode waves 



(a = 0.1, b = 2.0, c = ±2.5 in eq.(3.7), in a background spacetime with A'o = 0.025. We plot the L2 norm of the Hamiltonian 
and momentum constraints, C^^ and C^/^, for each cases. We see from all of them that strongly hyperbolic system improves 
the violation of the constraints. 
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plus-mode wave propagation 



- weakly hyp. 
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FIG. 5. Comparisons of the symmetric hyperbolic system (Ashtekar III) with the weakly hyperbolic system (Ashtekar 
original) (A^ = 1 slice), in the same way with Fig.^. We applied the same parameters with those of Fig.^. Figs, (al) and (a2) 
are of +-mode waves propagation and L2 norm of C^^^ and C^/^, respectively. Figs, (bl) and (b2) are of x-mode waves 

.ASH , ^ASH 



propagation and L2 norm of and Cm 

the violation of the constraints. 



respectively. We see from all of them that symmetric hyperbolic system improves 
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FIG. 6. Comparisons of 'adjusted' syste m wi th the different multiplier, k, in eqs. (5.1) and (5.2). The model was +-mode 



pulse waves (a — 0.1, b — 2.0, c — ±2.5 in eq.(3.6) in a background Kq = —0.025. Plots are of L2 norm of the Hamiltonian and 
momentum constraint equations, C^^^ and Cm [Figs-(a) and (b)] respectively. We see some ft produce better performance 
than the symmetric hyperbolic system. 
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FIG. 7. Comparisons of 'adjusted' system with the different multiplier, k, in eq. (B3). Plots are of L2 norm of the constraint 
equations, (Bl). 
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